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Abstract. Caesium fountain frequency-standards realize the second in the 
International System of Units with a relative uncertainty approaching 10~ 16 . 
Among the main contributions to the accuracy budget, cold collisions play 
an important role because of the atomic density shift of the reference atomic 
transition. This paper describes an application of the Bayesian analysis of 
the clock frequency to estimate the density shift and describes how the Baycs 
theorem allows the a priori knowledge of the sign of the collisional coefficient 
to be rigourously embedded into the analysis. As an application, data from the 
INRIM caesium fountain are used and the Bayesian and orthodox analyses are 
compared. The Bayes theorem allows the orthodox uncertainty to be reduced by 
28% and demonstrates to be an important tool in primary frequency-metrology. 
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1. Introduction 

Atomic fountains depend on laser cooling of caesium atoms down to a temperature 
of 1 fiK or an even lower value; this, together with the use of a single microwave 
cavity to implement the Ramsey separated-field spectroscopy, allowed the accuracy of 
primary frequency standards to be improved by more than one order of magnitude. 
Today, fountains realize the second with a relative accuracy ranging from 3 x 10 -16 
to 10 x 1CT 16 [HHIHHGEHEIIE] and ' in the last ten y ears > allowed the uncertainty 
of the Atomic International Timescale (TAI) unit to be reduced from 1 x 1CP 14 to 
5 x 1CT 16 . 

In Cs fountain clocks, cold collisions occur between the ultra-cold atoms with 
consequent perturbation of the atomic energy levels and shift of the atomic reference 
frequency. This shift, is proportional to the atom density and it must be carefully 
evaluated to correct the clock frequency. Since it is a major component of the 
uncertainty budget, it has been the subject of many theoretical and experimental 
studies [HEII. It depends on the collisional dynamics during the atom ballistic flight, 
which is related to the way the atom cloud is prepared before launch. When an atom 
cloud, captured into a magncto-optical-trap, is launched, its high initial density limits 
the collision energy and the density shift is strongly dependent on the temperature 
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and the mixture ratio [10] . However, if, as in fountains operation, the atom cloud is 
prepared with direct molasses capture or molasses expansion after magncto-optical- 
trap capture |T] , the collision energy is high during the whole flight and the shift has 
a negligible sensitivity to the atom temperature. Additionally, fountains are operated 
in a way that the collision energy is between 7 [iK and 10 fiK range and that the 
atomic sample is almost an evenly quantum mixture of the hyperfine eigenstates 
\F = 3,itif = 0) and \F = 4,mF = 0). In this case, the density shift has a negligible 
dependence on the mixture ratio and it is linearly dependent on density through a 
negative coefficient [3]. 

In a previous paper [11], we used the Bayesian inference [12l [13] to estimate 
the collisional coefficient, consistently with the constraint of a negative value. This 
estimate can be used to extrapolate the fountain frequency to the zero-density value. 
In the present paper we illustrate how to perform the extrapolation, given two or 
more frequency measurements at different atomic densities, irrespectively of the value 
of the collisional coefficient value, but consistently with its negative sign. From a 
mathematical viewpoint, the problem is to find the best-fit line through two or more 
points, with the constraint that the regression coefficient is negative. 

Our goal is twofold. First, we want to illustrate a non-trivial application of the 
Bayes theorem and the relevant data analysis. Second, we want to assess the viability 
a density shift estimation consistent with the constrain of a negative collisional 
coefficient and to test its performances. A Bayesian approach is important because 
the shift has the same magnitude as the measurement noise. Therefore, despite the 
fact that the shift is a linear function of the atom density having a negative regression 
coefficient, measurement results having the low-density frequency lower than the high- 
density one are relatively common. A linear extrapolation to zero density from these 
data is clearly meaningless. On the contrary, a Bayesian inference makes it possible to 
deal with this situation, thus avoiding physical absurdity and, consequently, reducing 
the extrapolation uncertainty. 

After the problem statement, in section 13.11 the Bayesian solution is given, 
first for two frequency measurements and, then, for list of frequency measurements. 
Eventually, in section l3~5l the Bayesian analysis is applied to data collected in a test 
measurement. All the symbolic and numerical calculations have been performed with 
the aid of Mathematica [14] . 

2. Experimental techniques 

The density shift is commonly evaluated by means of a differential measurement 
approach. The fountain is operated alternating low with high atom density and 
the frequency of a hydrogen-maser is measured using the fountain in the two 
configurations. By virtue of its frequency stability, in both short and medium time 
scales, the maser is used as a flywheel oscillator: the comparison of the two measured 
frequencies cancels the maser frequency and allows density shift to be evaluated. As 
reported in [9] , if the ratio between the atomic density and the total number of detected 
atoms is assumed constant, we can state that the density shift is proportional to 
the number of detected atoms. The differential measurement provides a collisional 
coefficient which can be used to extrapolate the frequency data to zero density. 

The actual experimental practice in the various laboratories differs in the way 
the atom clouds of different densities are prepared and in the durations of the 
fountain operation at low and high density. The extrapolations techniques of the 
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fountains of the Physikalish Tcchnischc Bundesanstalt (PTB-CsFl), of the National 
Physical Laboratory (NPL-CsFl), and of the National Institute of Standards and 
Technology (NIST-F1) differ only by the driving parameter and the time scheduling 
of the differential measurements. At the Sytemcs de Reference Temps Espacc, a rapid 
adiabatic passage technique is used, which ensures very accurate evaluation of the 
collisional coefficient because the atom density is precisely set at one value and at its 
half. 

At INPJM, we vary the atom density through the loading time of the magneto- 
optical trap. To evaluate the density shift, the fountain is operated alternately at 
high and low atomic density; the 70 ms and 300 ms loading times provide respectively 
the low and high density configurations. The ratio between the number of detected 
atoms in the two configurations ranges from three to four. Since in the low density 
configuration stability is poor and resolution is limited by the atom shot noise, 
the fountain is operated alternating about 21000 s in the low density configuration 
and about 6000 s in the high density one. The hydrogen maser frequency is then 
extrapolated to zero density by 



where y is the sought zero density frequency, y\ and y<i are the frequency in low and 
high density conditions, and x\ and x-i are the number of atoms detected at the low 
and high density, respectively. The extrapolation is carried out for each (yi, j/2) pairs. 
The total duration of each run is 27000 s; this ensure a good rejection of systematic 
effects - fluctuations of the hydrogen maser frequency, magneto-optical trap, and atom 
detection efficiency - which could bias the frequency extrapolation. 

By neglecting the a x uncertainty of the atom-density measurements - an 
acceptable omission provided a x <C <J y j\yi — y\\ - the uncertainty of ((T|) is 



where a y is the uncertainty of the frequency measurements. Deviations from a 
linear relationship and an imperfect rejection of long term effects are estimated to 
contribute to the uncertainty by 20% of the density-shift value. Here, this non- 
statistical contributions to the uncertainty have not been considered, though they take 
part of the total uncertainty budget of atomic clocks. Besides, being the statistical 
contribution of the same order of magnitude as the shift, model errors do not contribute 
significantly. 

3. Zero density extrapolation 

3.1. Statement of the problem 

In the simplest case, we want to perform a linear regression analysis of two clock- 
frequency, at low and high atom densities, given the prior knowledge that the slope of 
the sought line is negative. In particular, we are interested in the value of the intercept 
of the regression line. The data arc assumed normally distributed about axi + b with 
the same standard deviations a yi that is, 



a = 
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The joint probability density of the data is 

rw I i\ 1 r (yi - axi - b) 2 + (y 2 - ax 2 - b) 2 

Pv(Vi,y2\a,b) = ^2 ex P L ^2 J- ( 4 ) 

The question is how to account for the a < constraint and to infer a consistent value 
for b. 

Here and in the following, we will use the notation P r (r,|sj) to indicate the 
probability density that the quantity r has the particular value r = ri if the 
parameter s in the probability distribution has the particular value s = Sj. Irrelevant 
conditionals, such as Xi and a y in (J3J), will be dropped. 



3.2. Bayesian solution 

The Bayesian approach takes account of a < into the probability density -P a fc(ao> bo) 
of the regression coefficient and intercept values before the measurement result is 
know; in the absence of any additional information let us assume that 

P a b{a ,b ) = #(-oo)) (5) 

where $(.) is the Heavisidc function. According to the Bayes theorem the post-data 
probability density of the a and b values is 

P a b(ao,bo\yi,y2) oc Pu(yi,y2\ao,b )P a b{ao,bo) 

r (yi - aoxi - bo) 2 + (j/2 - a x 2 - bo) 2 , 

ocexp|_ — Jtf(-ao), ( 6 ) 

za y 

where a meaningless normalization coefficient has been omitted. This probability 
density embeds both the a < constraint and the information delivered by the data. 

To obtain the probability density of the zero-density frequency, irrespectively of 
the regression coefficient, oq is integrated out from ([6]) by marginalization. Hence |14j . 



Pb{bo\yi,V2) 
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Figure 1. Probability density of the intercept values; 212/2:1 = 3, solid (red) 
is (3/2 — Vl)/ a y = — 4, dashed (blue) is (j/2 — yi)/<?y = 0, dotted (black) is 
(s/2 - yij/o-y = +4. 
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where y and a 2 are the intercept of the line through the data and its variance - given 
by HD and ©, 

y = , (8) 

x\ + x 2 

is the mean of the measured frequencies weighed by the relevant atomic densities, 

2 x\ + X 2 2 (a \ 
a =7 so °u ( y J 

is the uncertainty of y, and erfc(.) is the complementary error function. The post- 
data probability density is shown in Fig. [TJ The Gaussian factor is the probability 
density we obtain by a classical analysis: it is the probability density of the orthodox 
extrapolation, having y mean and a 2 variance. The erfc factor, originating from the 
application of the Bayes theorem and from marginalization, takes account of the a < 
constraint and cuts off the intercept values less than y. 

The probability density {7} is the result of the Baycsian analysis. To convert 
it into a single numerical estimate, a loss associated with the estimate error must 
be specified. The optimal estimator minimizes the expected loss over (JT)). A loss 
proportional to the squared or absolute errors indicates the mean or the median, 
respectively; a constant loss indicates the most probable values. Confidence intervals 
are easily expressed by integrating {TJ to obtain the cumulative distribution function. 

Figure [5] shows the mean and standard deviation of the zero-density frequency. 
In the limit case when {y 2 — V\)l a y ^ 0, the mean tends to the y intercept of the line 
through the data, but with larger uncertainty. This is explained by observing that, 
via marginalization, the Bayesian analysis takes account of all the possible regression 
coefficients. The understanding of the [y 2 — yi)/<^ y S> limit requires some training 
of our intuition. Let us observe that 14 



lim P b (b \yi,y 2 ) = -j= exp 
{V2-yi)/^y^oo V 7ra y 



(bo-m (10) 



where y = (yi + y 2 )/2. Therefore, the intercept mean approaches the sample mean. 
In this case, the data are inconsistent with the a priori knowledge on a, but we did not 
allow room for such inconsistency when formulating the problem and, though at least 
one of the data is clearly wrong, we do not know which it is. Consequently, the best we 
can do is to fit the data with a line satisfying the a < constraint. The sample mean 
corresponds to a = 0; the Bayesian inference is slightly greater because it accounts 
for all the a < values. When (y 2 — yi)/<J y 3> 0, the extrapolation uncertainty 
is smaller than the uncertainty of the extrapolation when (jj 2 — yi)/&y *C 0. The 
reason is that, when (y 2 — yi)/cr y — > 00, we did not question the <j y /\/2 uncertainty 
of the sample-mean, regardless of data inconsistency. Finally, we observe that, when 
(j/2 — V\)l a y ~ 0, the mean is always greater than y\ and smoothly connects the y 
and y asymptotes. 

When (2/2 — yi)/(j y w 0, the Bayesian analysis seems to overestimate the zero- 
density frequency, which, from an orthodox analysis, is expected to be quite near 
2/1. The supposed overestimation is due to the scarce information delivered by the 
data and the use of ([S]): with a uniformly negative regression coefficient, a vertical 
regression is, a priori, as probable as a horizontal one and the Bayes theorem accounts 
for both possibilities. In the appendix A, we give the result of a Monte Carlo simulation 
which confirms that ^ describes correctly the pre-data probability density, when no 
information, apart a < 0, is available. To summarize the results, we used the mean 
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Figure 2. Mean and standard deviation of the zero-density frequency, given 
{xi = 1, j/i } and {12 = 3, 1/2}- The left and right sides of the diagram correspond 
to physical and unphysical values of the collisional coefficient, respectively. 
The coloured area indicates the extrapolation uncertainty. The straight lines 
(asymptotic limits) are b = y\ — (1/2 — yi )/2 (best-fit line intercept) and b = 
(yi + j/2)/2 (sample mean). Dashed line is the mean when the same {l,yi} and 
{3,1/2} pairs are observed three times. 



and the standard deviation as calculated from (JT]). However, when (j/2 — 2/i)/ (J y ~ 0, 
the post-data probability density is distorted; in this case, the most probable value is 
somewhat smaller than our estimate. 

The role of the assumed prc-data probability density of the regression coefficient 
is also shown by considering the Bayesian estimate of the zero-density frequency given 
a single data pair, for instance, {x\ = 1, y\ = 0}. In this case, the marginal probability 
density of the best-fit line intercept is the improper distribution 

P 6 (&o|zi)aerfc(--^-J, (11) 

where roughly all the frequencies greater than the measured value arc equally probable. 
However, as an increasing number of data become available, the post-data distribution 
is dominated by the likelihood function; P a b(ao, bo) becomes irrelevant, and we are led 
to the same conclusion irrespectively of the prior probability density. 



3.3. Extension to repeated observations 



To extend the previous analysis to N measurement pairs {xi, yi}, we must rewrite the 
joint probability density of the data as 



N 



P»{y\a, b) = Y[ N(y t \ax t + b, a 4 2 ), 



(12) 



where y is the vector of the measured frequencies, N(yi\i/i, erf) is a normal distribution 
with Vi mean and af variance, and we assumed the data independent and identically 
distributed. Hence, the joint post-data probability density of the regression coefficient 
and zero-density frequency is 

Pab(ao,b \y) oc P v (y\a, 6)i?(-a ) 



oc exp 



( a - a b - y ) C ab 



ap — a 

bo-y 



(13) 
0(-ao), 
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where, as it is shown in the appendix B, a, y, and C a b are the least-squares estimates 
and covariance matrix of a and b. Eventually, in the same way as in , the parameter 
ao is integrated out of the problem by marginalization. There is no additional insight 
in integrating (|13p analytically, but it is must be noted that, when a/a a <C 0, where 
a a is the least-squares uncertainty of a, then E(b) — > y, where E(b) is the b mean over 
the Pb(bo\xi, y%) distribution. On the contrary, when a/a a 3> 0, then E(b) — * y, where 
y is the sample average of the data. In the general case, the smaller the uncertainty 
of the least-squares line is, the closer is the extrapolated frequency E(b) to y, if a < 0, 
and to y, if a > 0. This is shown in Fig. [5] for the particular case when the same 
and {3,j/2} pairs are observed three times. 

S.^. Discussion of results 

We can also process the data pairs sequentially In this case, we start with ([6]) based 
on the first data pair. This probability density substitutes for ([5]) as the pre-data 
probability density in the analysis of the second data pair. When this procedure 
is repeated up to the last pair, we obtain the same post-data distribution as that 
obtained with the one-step approach. To better understand the Bayesian analysis, let 
us consider three identical data-pairs, where (2/2 — y\)/°~y is close to zero. By starting 
with ([S]), given the first data pair, and by repeating the analysis given the second 
and, then, given the third, we note that the Bayesian extrapolation changes each time 
a new pair is made available. This result is in striking contrast with the orthodox 
analysis where, if we use three identical data pairs, the extrapolation uncertainty is 
reduced, but the extrapolated value remains the same. 

This apparent paradox is explained by observing that, when applying the Baycs 
theorem for the first time, the only piece of information is that the regression 
coefficient value is between zero and minus infinity with equal probability. Hence, 
the extrapolation is biased towards a frequency value much greater than y\. In the 
second iteration, additional information is available, namely the result of the first 
measurement. Both pieces of information are synthesized in the post-data probability 
density ((6j, which substitutes for the pre-data probability density ([5|) and limits 
extrapolation to a neighbour of the classical one. In each subsequent application of 
the Bayes theorem, we update the pre-data probability density, thus further reducing 
the discrepancy between orthodox and Bayesian extrapolations. The Bayes theorem, 
by prescribing that the regression coefficient is negative and that all its negative 
values including those arbitrarily large - have the same prior probability, infers 
that the classical extrapolation statistically underestimates the zero-density frequency. 
However, when several measurement results are available, and all are consistent with 
a < 0, Bayesian extrapolation approaches to the classical one. 

In the absence of informative data, the sensitivity of the Bayes theorem to prior 
information, synthesized in the pre-data probability density of the measurand, may 
be disappointing. This could discourage the use of Bayesian methods, in view of an 
apparent lack of objectivity. However, a seminal paper by R. T. Cox |15| demonstrates 
that, in order to make consistent inferences, it is necessary to resort the Bayes theorem. 
If we give it up, because we are adverse to make an estimate depending on the prior 
information, we put our results at the risk of contradictions. 
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3. 5. Estimation of the hydrogen maser drift 

Given its good stability during time intervals from hours to weeks, the hydrogen maser 
is a common choice as a frequency flywheel for the differential measurements. It is also 
used as a transfer oscillator when the fountain is used to evaluate the TAI time unit. 
However, the hydrogen-maser frequency drifts linearly; this drift is usually evaluated 
and removed by means of orthodox statistical techniques. Since the drift is stable for 
long time intervals, we can use the knowledge accumulated in the past fountain run to 
exploit a full Baycsian simultaneous evaluation of the drift and zero-density frequency, 
with the use of a two-dimensional linear regression. 

This requires a slight modification of the sampling distribution of the data, which 
are therefore assumed normally distributed about ax, + b + cfcj, where Xi and ij are 
respectively the density and epoch related to the frequency value yi , a is the collisional 
coefficient, c is the hydrogen maser drift and b is the frequency value extrapolated to 
zero density and the epoch t = 0. Hence, the joint probability density of the data 
becomes 

N 

P v (y\a, b, c) = Yl N( yi \ax t + b + ct z , aj). (14) 

i=l 

The drift of the hydrogen maser frequency is known; therefore, its pre-data probability 
density is 

P c (c ) =N(co\c,a c ), (15) 

where c ± o~ c is its pre-data estimate. Eventually, the joint post-data probability 
density of the model parameters is 

N 

Pabc(ao,b ,c \y) cx _Q N(yi\a Xi + b + c U, af)N(c \c, cr c )i?(-a ) (16) 

?:=i 

and the post-data probability densities of each individual parameter irrespectively of 
the others are obtained by marginalization. Then, the marginal post-data probability 
density of the zero-density frequency is 

/+oo 
Pabc{ao,b , co\y) da dc (17) 
-oo 



4. Frequency extrapolation in a Cs fountain 

Tabic I records measurement results and the relevant uncertainties, collected during a 
TAI unit evaluation run of the IT-CsFl fountain at INRIM. The density values have 
been so scaled that the mean low-density is unitary; the frequency values have been 
given in units of the a v = 3.9 x lCP 15 z^cs uncertainty - where z^cs = 9 192 631 770 Hz 
- and have been so shifted that the mean low-frequency is one. The same data are 
shown in Fig. [3] 

Some of the low- and high-density data pairs in Table 1 lie on a positive-slope 
line, in agreement with Gaussian dispersion of the frequency data. By the orthodox 
approach, the prior information is not taken into account and these positively-sloped 
data contribute to the best-fit line with the same weight as the others. It is not so 
when a Bayesian analysis is made. By applying the model described in !3.5l to the data 
in Table 1, the marginal post-data probability density of the zero-density frequency 
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Table 1. Clock frequency vs. atom density. Measurement have been performed 
alternately at low and high densities at equal 0.313 day intervals; time increases 
from top to bottom and, then, from left to right. p low and i/ low are the mean low- 
density and -frequency, a v = 3.9 X 10 i^Cs is the mean uncertainty of frequency 
data. 
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Figure 3. Results of low (red squares) and high (blue bullets) density frequency 
measurements. The hydrogen-mascr drift has been removed and the line is the 
intersection of the best-fit plane with t = 0. 



is given by (fl7]) and it is shown in Fig. [U where we used the c = (0.41 ± 0.05)cr„/day 
pre-data estimate of the hydrogen mascr drift. 

After the relevant marginalization, the mean and standard deviation of the 
hydrogen maser drift, collisional coefficient, and zero-density frequency have been 
calculated and are shown in Table 2 [T3]. Since, in general, the probability densities 
of these quantities are not Gaussian, the mathematical definitions of the mean and 
variance, for instance, 

/+oo 
d£ (18) 
-oo 

and 

var(6)= / [£-E(&)] 2 P 6 (£|y)d£, (19) 

«/ — OO 
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have been used. The same table shows also the classical estimates obtained by fitting 
the data with the axi + b + ctj plane. 

When comparing the Baycsian and the orthodox frequency extrapolation, we 
observe that the Bayes theorem allows the orthodox uncertainty to be reduced by 28%. 
This is in agreement with the best usage of the prior information, which usage reduces 
significantly the variability of the regression coefficient. In terms of original relative 
frequency units, the Baycsian analysis allows the extrapolation uncertainty to be 
reduced from 9.8 x 10 _16 i/cs to 7.0 x 10 _16 ^cs- We also checked that the use of a more 
conservative pre-data estimate of the hydrogen maser drift, e.g., (0.50 ± 0.15)cr„/day, 
does not significantly affect the post-data probability density. 



Table 2. Comparison of classical and Baycsian estimates of the hydrogen maser 
drift, collisional coefficient, and zero-density frequency. a v = 3.9 X 10 _15 i^c s is 
the mean uncertainty of frequency data 





classical analysis 


Bayesian analysis 


H maser drift 
collisional coefficient 
zero-density frequency 


(0.43 ± 0.03)ov/day 
(-0.14±0.10)ov/p low 
z/ low + (0.18 ±0.25)cr„ 


(0.44 ± 0.02)ov/day 
{-0.18±0m)a u /p low 
^ 1ow + (0.24±0.18)ct„ 



A final consideration concerns the conceptual interpretation of classical and 
Bayesian results. When summarized by the orthodox best estimate and uncertainty, 
the available information about the zero-density frequency is synthesized by the 
realization of a random variable (the extrapolated frequency value) and by a measure 
of the width of distribution from which this realization has been sorted - e.g., the 
variance. On the contrary, the Bayesian analysis synthesizes the information by a 
measure of the measurand-distribution location and by a coverage interval containing 
the value of the zero-density frequency with stated probability. 

5. Conclusions 

When extrapolating the frequency of a caesium-fountain clock to zero atom-density, 
the Bayes theorem makes it possible to take account of a negative correlation between 



2.5 




-1.0 -0.5 0.0 0.5 1.0 

(V0 - Vlow)0-y 



Figure 4. Post-data probability density of the zero-density extrapolated 
frequency. 
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the density and the frequency at the very beginning of data analysis, in a way ensuring 
logical consistency of the extrapolation. 

In the specific example here considered, the theorem demonstrates itself capable 
of estimating the hydrogen maser drift and the collisional coefficient at the same time 
and of improving the orthodox extrapolation of the zero density frequency With use of 
the sign information, the frequency extrapolation can be efficiently carried out also in 
the case of a single data pair and also when the measured frequencies are inconsistent, 
that is, when they violate the yi > y\ constraint. 

We have considered the accuracy evaluation of the INRIM Cs fountain during 
standard operation; the uncertainty contribution due to atomic density shift is reduced 
by 28%, from 9.8 x 10~ 16 to 7.0 x 10~ 16 . In this application, the frequencies are 
averaged classically over 21000 s (low-frequency configuration) and 6000 s (high- 
frequency configuration); the Bayesian analysis is then applied to the resulting set 
of data. Since these averages do not account for a negative collisional coefficient, 
future work will be aimed at investigating the optimal split the 27000 s block into 
shorter ones to reduce the influence of the classical averaging. 

Appendix A. Inverse Monte Carlo simulation 

The probability density of the zero-density frequency f7]) was calculated also 
numerically by inverse Monte Carlo simulation. In direct Monte Carlo simulation 
a random list of measurement results, having a fixed measurand value, is generated 
by repetitions of a numerical experiment. On the contrary, in inverse simulation it is 
generated a random list of measurand values having a fixed measurement result. 

A brute force approach exemplifies the simulation procedure. Let the measured 
values of the low- and high-density frequencies be fixed, for example, y\ = and 
i/2 = — 1 ~ in the normalized units where <j y = 1. In the inverse simulation, y\ and y<i 
are sampled from N(i/i\axi + b, 1), where the slope and intercept are randomly chosen 
according ([5|). If the sampled frequencies are y\ = and j/2 = —1, that is, if they 
happen to be identical to the wanted measurement results, the sorted zero frequency 
b is appended to the Monte Carlo list. 

This procedure is thoroughly inefficient. A short-cut is to observe that any 
arbitrary {y±, j/ 2 } pair - sorted from N(jji\aXi + b, 1), with arbitrarily chosen values of 



0.4 




O'o-;vi)/o\v 



Figure 5. Monte Carlo histogram of the zero-density frequency, given {x± = 
1,2/1 = 0} and {x2 = 3, 3/2 = — !}■ Solid line is the theoretical prediction J7J| 
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a and b - can be mapped into the wanted {0, —1} pair by 

(1 + 2/2 -yi)(xi -xi) 
Vi^Vi yx- A.i) 

x 2 - xi 

The same result should have been obtained if the data pair should have been sampled 
from N(i>i\vi, 1), where 

/ 1 + 2/2 - 2/i \ , , . 1 + Vi - 2/1 , . „s 

i/j = a + o H X\ — y\. 

V xi- x x J x 2 - xi 

Hence, given any {yi-,y 2 } sampled from N{yi\axi + 6,1), the inverse Monte Carlo 
simulation sorts randomly the zero-frequency 

, , . 2/2 - yi Q s 
^ N xi-yi, (A.3) 

which is obtained by setting Xi — in (|A.2[) . Provided the relevant collisional 
coefficient in (|A.2|) . a — (1 + y 2 — yi)/(x2 — a^i), is negative - as requested, (|A.3[) 
is appended to the Monte Carlo list, otherwise it is rejected. 

A Mathcmatica script illustrating the inverse Monte Carlo simulation is appended 
below; the simulation results are shown in Fig. the Monte Carlo frequencies agree 
with the marginal probability density of the zero-density frequency predicted by ([7]). 
We are unable to find out if and where the uniformity of the pre-data distribution has 
been used in the simulation; therefore, this result confirms that ([5]) indicates correctly 
the absence of any prior information, apart from the sign. 

(* Inverse Monte Carlo Simulation *) 

np = 100000; (* # Monte Carlo runs *) 

a = -0.5; b = 0; (* model parameters, initially -0.5 and *) 
xl = 1; x2 = 3; (* low and high densities *) 
(* list of low density frequencies *) 

yl = a xl + b + RandomReal [NormalDistribution [0 , 1], {np}] ; 
(* list of high density frequencies *) 

y2 = a x2 + b + RandomReal [NormalDistribution [0 , 1], {np}] ; 

(* list of updated collisional coefficients *) 

a = a - (1 + y2 - yl)/(x2 - xl) ; 

(* list of updated zero frequencies *) 

b = b + (1 + y2 - yl) xl/(x2 - xl) - yl; 

ab = Transpose [{a, b}] ; 

ab = Select[ab, #[[1]] < &] ; (* selects a < values *) 
b = Transpose [ab] [[2]] ; (* Monte Carlo list *) 

Appendix B. Sufficient statistics for the zero-density frequency and 
collisional coefficient 

The post-data probability density (fT3")) holds because a and y - the slope and 
intercept of the best-fit line through the data - are sufficient statistics for the 
collisional coefficient and zero-density frequency. To demonstrate this, let the sampling 
distribution (fT2l) be written as 



P„(y|/3) oc exp 



h y -Afl) T C-\ V -A0) 



(B.l) 
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where y is the vector of the measured frequencies, (3 is the vector of the unknowns, A 
is the design matrix, and C y the covariance matrix. The exponent of (|B.2[) . x 2 (/3), is 
a quadratic form in /3; hence 



where p = CpA T C~ 1 y is the least-squares estimate of (3 and Cp = (A T C y ~ 1 A)~ 1 its 
covariance matrix. Eventually, by leaving out the terms independent of /3, which are 
unessential in (fT3|) . 
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P v (y\0)ocexp --(p-pfC^(p-p) 
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